﻿#ifndef TURBO_MODEL
#define TURBO_MODEL

#include <map>
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include "grid.h"
#include "datatypes.h"
#include "optimization.h"
#include "interpolation.h"


//Step info
class StepInfo {
public:
	double TimeStep;
	double Residual;
};


//Medium properties
class MediumProperties {
public:
	double Gamma;
	double Cv;
	double Viscosity;
};




//Class that represents model and its underlying properties
template<class RiemannSolver, class TurbulentModel> 
class Model {
public:			
	//Boundary conditions section
	//Boundary condition types availible in model
	enum BCType{
		SubsonicInlet,
		SubsonicOutlet,
		NoSlip,
		Symmetry,
		Farfield
	};		

	class BoundaryCondition {		
		Dictionary parameters;
	public:
		Model& model;
		BCType type;
		
		void setParameter(std::string name, double value) {
			parameters.setValue(name, value);
		};	

		BoundaryCondition(Model& _model) : model(_model) {}; 

		virtual ConservativeVariables getDummyValues(ConservativeVariables UL, const Face& face) = 0;
	};


	//Boundary conditions availible
	class NoSlipBoundaryCondition : public BoundaryCondition {
	public:
		NoSlipBoundaryCondition(Model& _model) : BoundaryCondition(_model) {};

		ConservativeVariables getDummyValues(ConservativeVariables inV, const Face& face) {
			inV.rou *= -1;
			inV.rov *= -1;
			inV.row *= -1;
			return inV;
		};
	};

	class SymmetryBoundaryCondition : public BoundaryCondition {
	public:
		SymmetryBoundaryCondition(Model& _model) : BoundaryCondition(_model) {};

		ConservativeVariables getDummyValues(ConservativeVariables inV, const Face& face) {
			Vector roU(inV.rou, inV.rov, inV.row);			
			Vector newRoU = roU - 2*(roU * face.FaceNormal)*face.FaceNormal/face.FaceNormal.mod();
			inV.rou = newRoU.x;
			inV.rov = newRoU.y;
			inV.row = newRoU.z;
			return inV;
		};
	};

	class SubsonicInletBoundaryCondition : public BoundaryCondition {
		double _pressure;
		double _temperature;
		Vector _velocity;
		MediumProperties medium;
	public:		
		SubsonicInletBoundaryCondition(Model& _model) : BoundaryCondition(_model) {			
			medium = model.medium;
		};

		void setParams(double pressure, double temperature, Vector velocity) {
			_pressure = pressure;
			_temperature = temperature;
			_velocity = velocity;
		};				

		ConservativeVariables getDummyValues(ConservativeVariables inV, const Face& face) {			
			//Compute outgoing riemann invariant
			//Obtain speed of sound			
			Vector vd = model.GetVelocity(inV);
			double cd = model.GetSoundSpeed(inV);
			double Rminus = vd * face.FaceNormal / face.FaceNormal.mod() -  2*cd/(medium.Gamma - 1.0);
			//Now solve equation Rminus = const for velocity
			//Use Blazek recomendations	
			/*double cost = -1.0;
			if (vd.mod() > std::numeric_limits<double>::epsilon()) cost = - vd * face.FaceNormal / (face.FaceNormal.mod() * vd.mod());
			double cd2 = cd*cd;
			double c02 = cd2 + (medium.Gamma - 1.0) * vd * vd / 2.0;
			double A = (medium.Gamma - 1.0) * cost * cost + 2;
			double cb = A * c02 / (Rminus * Rminus * (medium.Gamma - 1.0)) - (medium.Gamma - 1.0) / 2.0;
			cb = sqrt(cb);
			cb = 1 + cost * cb;
			cb *= (-Rminus *  (medium.Gamma - 1.0)) / A;
			double cb2 = cb * cb;
			double T0 = _temperature;	
			double k = cb2/c02;
			if (k > 1.0) k = 1.0;
			double Tb = T0 * k;
			double P0 = _pressure;
			double Pb = P0 * pow(Tb/T0, medium.Gamma / (medium.Gamma - 1.0) );
			double vbmod = sqrt(2* medium.Gamma * medium.Cv * (T0 - Tb));
			Vector vb = _velocity * vbmod / _velocity.mod();*/		
			Vector vb = _velocity;
			double P0 = _pressure;
			double Pb = P0;	
			double cb = (vb * face.FaceNormal / face.FaceNormal.mod() - Rminus);
			cb = cb * 0.5 * (medium.Gamma - 1.0);
			double eb = cb*cb / ((medium.Gamma - 1.0) * medium.Gamma);
			double Tb = eb / medium.Cv;
			return model.PrimitiveToConservativeVariables(vb, Pb, Tb, medium);
		};
	};

	class SupersonicInletBoundaryCondition : public BoundaryCondition {
		double _pressure;
		double _temperature;
		Vector _velocity;
		MediumProperties medium;
	public:		
		SupersonicInletBoundaryCondition(Model& _model) : BoundaryCondition(_model) {			
			medium = model.medium;
		};

		void setParams(double pressure, double temperature, Vector velocity) {
			_pressure = pressure;
			_temperature = temperature;
			_velocity = velocity;
		};				

		ConservativeVariables getDummyValues(ConservativeVariables inV, const Face& face) {												
			double Pb = _pressure;
			double Tb = _temperature;			
			Vector vb = _velocity;			
			return model.PrimitiveToConservativeVariables(vb, Pb, Tb, medium);
		};
	};

	class InletBoundaryCondition : public BoundaryCondition {
		double _pressure;
		double _temperature;
		Vector _velocity;
		MediumProperties medium;

		//Boundary condition to choose
		SupersonicInletBoundaryCondition supersonicBC;
		SubsonicInletBoundaryCondition subsonicBC;
	public:		
		InletBoundaryCondition(Model& _model) : BoundaryCondition(_model), supersonicBC(_model), subsonicBC(_model) {			
			medium = model.medium;			
		};

		void setParams(double pressure, double temperature, Vector velocity) {
			_pressure = pressure;
			_temperature = temperature;
			_velocity = velocity;
			supersonicBC.setParams(pressure, temperature, velocity);
			subsonicBC.setParams(pressure, temperature, velocity);
		};				

		ConservativeVariables getDummyValues(ConservativeVariables inV, const Face& face) {												
			double cd = model.GetSoundSpeed(inV);
			Vector vd = model.GetVelocity(inV);
			if (cd < -vd * face.FaceNormal / face.FaceNormal.mod()) {
				return supersonicBC.getDummyValues(inV, face);
			} else {
				return subsonicBC.getDummyValues(inV, face);
			};
		};
	};


	class SubsonicOutletBoundaryCondition : public BoundaryCondition {
		double _pressure;
		MediumProperties medium;
	public:
		SubsonicOutletBoundaryCondition(Model& _model) : BoundaryCondition(_model) {			
			medium = model.medium;
		};

		void setParams(double pressure) {
			_pressure = pressure;
		};		

		ConservativeVariables getDummyValues(ConservativeVariables inV, const Face& face) {												
			//Now solve equation Rminus = const for velocity
			//Use Blazek recomendations						
			double Tb = model.GetTemperature(inV);			
			double Pb = _pressure;			
			Vector vb = model.GetVelocity(inV);			
			return model.PrimitiveToConservativeVariables(vb, Pb, Tb, medium);
		};
	};

private:
	//Internal variables
	Grid _grid;	//Grid
	std::map<int, BoundaryCondition*> _boundaryConditions; //Boundary conditions	

	//Store conservative variables for each cell
	DistributedEntityManager<ConservativeVariables> U;	

	//Riemann solver object
	RiemannSolver rSolver;

	//Tubulent model object
	TurbulentModel turbulentModel;

	//Wall information stored in each cell
	struct CellWallInfo {
		int wallFaceIndex;
		double angle;
		double distance;
	};

	//Information stored in each wall face
	struct FaceWallInfo {
		int wallFaceIndex;
		std::vector<int> layerCells;
		double boundaryLayerThickness;
	};

	alglib::kdtree kdtWall; 	//KDTree structure and related operations for wall faces
	alglib::kdtree kdtCells; 	//KDTree structure and related operations for all cells

	std::map<int, CellWallInfo> _wallInfo;	//cellIndex -> wall information
	std::map<int, FaceWallInfo> _wallFaces; //wall faceIndex -> wall information
	std::set<int> _wallBoundaryMarkers;

	//Computational settings
	double CFL;	//CFL condition value (0.35 for example)

	//Computed variables for each face
	std::map<int, std::vector<double>> fluxes;
	std::map<int, std::set<int>> cellNeigbours;
	std::map<int, std::vector<Vector>> gradsCells;	
	std::map<int, std::vector<Vector>> gradsFaces;	

public:	
	//Medium properties
	MediumProperties medium;

	//Properties
	StepInfo stepInfo;
	double totalTime;

	//Accessors and settings
	void SetGamma(double gamma) {
		medium.Gamma = gamma;
		rSolver.SetGamma(gamma);
	};

	void SetCv(double cv) {
		medium.Cv = cv;
	};

	void SetViscosity(double v) {
		medium.Viscosity = v;
	};

	void SetCFLNumber(double cfl) {
		CFL = cfl;
	};

	void SetHartenEps(double eps)
	{
		rSolver.SetHartenEps(eps);
	};

	//Constructor
	Model() {	
		totalTime = 0; //Init simulation time
	};

	void BindGrid(Grid& grid) {
		//Grid
		_grid = grid;

		//Allocate memory for data over each cell		
		std::vector<Cell*> cells = _grid.cells.getLocalNodes();
		for (int i = 0; i<cells.size(); i++) {					
			U.add(ConservativeVariables(cells[i]->GlobalIndex));
		};
	};

	void SetBoundaryCondition(std::string BCName, BoundaryCondition& bc) {
		if (_grid.patchesNames.find(BCName) == _grid.patchesNames.end()) throw Exception("No such boundary availible");
		int bcMarker = _grid.patchesNames[BCName];
		_boundaryConditions[bcMarker] = &bc; 
	};
	
	void SetInitialConditions(ConservativeVariables initValue) {
		std::vector<Cell*> cells = _grid.cells.getLocalNodes();
		for (int i = 0; i<cells.size(); i++) {
			U[cells[i]->GlobalIndex] = initValue;
			U[cells[i]->GlobalIndex].GlobalIndex = cells[i]->GlobalIndex;
		};
		return;
	};

	void SetWallBoundary(std::string BCName, bool isWall) {
		//Obtain boundary bc marker
		if (_grid.patchesNames.find(BCName) == _grid.patchesNames.end()) throw Exception("No such boundary availible");
		int bcMarker = _grid.patchesNames[BCName];

		if (isWall) {
			_wallBoundaryMarkers.insert(bcMarker);
		
			//Move through all faces and add wall face if needed
			for (std::set<int>::iterator it = _grid.patches[bcMarker].faces_idx.begin(); it != _grid.patches[bcMarker].faces_idx.end(); it++) {
				_wallFaces[*it] = FaceWallInfo();
				_wallFaces[*it].wallFaceIndex = *it;
			};
		} else {
			if (_wallBoundaryMarkers.find(bcMarker) != _wallBoundaryMarkers.end()) _wallBoundaryMarkers.erase(bcMarker);

			//Move through all faces and remove wall face if needed
			for (std::set<int>::iterator it = _grid.patches[bcMarker].faces_idx.begin(); it != _grid.patches[bcMarker].faces_idx.end(); it++) {
				if (_wallFaces.find(*it) != _wallFaces.end()) _wallFaces.erase(_wallFaces.find(*it));
			};
		};
	};

	void Step() {
		StepInfo info;
		info.TimeStep = std::numeric_limits<double>::max();					
		std::vector<Face*> faces = _grid.faces.getLocalNodes();		

		//Compute convective flux for each cell face and apply boundary conditions						
		fluxes.clear();
		for (int i = 0; i<faces.size(); i++) {
			Face& f = *faces[i];
			ConservativeVariables UL = U[f.FaceCell_1];
			ConservativeVariables UR;
			if (f.isExternal) {
				UR = GetDummyCellValues(UL, f);
			} else {
				UR = U[f.FaceCell_2];
			};
			std::vector<double> flux = rSolver.ComputeFlux( UL,  UR, f);

			//Store fluxes
			fluxes[f.GlobalIndex] = flux;

			//Adjust timestep
			double ts = _grid.cells[f.FaceCell_1].CellVolume/rSolver.MaxEigenvalue;
			if (info.TimeStep > ts) info.TimeStep = ts;
			if (!f.isExternal) {
				ts = _grid.cells[f.FaceCell_2].CellVolume/rSolver.MaxEigenvalue;
				if (info.TimeStep > ts) info.TimeStep = ts;
			}
		};

		//Compute gradients of conservative variables for each face
		ComputeGradients();

		//Compute viscous fluxes
		ComputeViscousFluxes();

		//Compute distances to the wall
		ComputeWallDistances();

		//Apply CFL condition
		info.TimeStep *= CFL;
		info.Residual = 0;

		//TO DO remove
		//Save previous values
		/*std::map<int, ConservativeVariables> deltaU;
		std::set<int> indxs = U.getAllIndexes();
		for (std::set<int>::iterator it = indxs.begin(); it != indxs.end(); it++) {
			for (int k = 0; k<5; k++) deltaU[*it][k] = 0;
		};*/


		//Distribute fluxes to cells
		for (int i = 0; i<faces.size(); i++) {
			Face f = *faces[i];
			std::vector<double>& flux = fluxes[f.GlobalIndex];						

			U[f.FaceCell_1].ro -= info.TimeStep * flux[0]/_grid.cells[f.FaceCell_1].CellVolume;
			U[f.FaceCell_1].rou -= info.TimeStep * flux[1]/_grid.cells[f.FaceCell_1].CellVolume;
			U[f.FaceCell_1].rov -= info.TimeStep * flux[2]/_grid.cells[f.FaceCell_1].CellVolume;
			U[f.FaceCell_1].row -= info.TimeStep * flux[3]/_grid.cells[f.FaceCell_1].CellVolume;
			U[f.FaceCell_1].roE -= info.TimeStep * flux[4]/_grid.cells[f.FaceCell_1].CellVolume;			

			//TO DO remove
			/*deltaU[f.FaceCell_1].ro -= info.TimeStep * flux[0]/_grid.cells[f.FaceCell_1].CellVolume;
			deltaU[f.FaceCell_1].rou -= info.TimeStep * flux[1]/_grid.cells[f.FaceCell_1].CellVolume;
			deltaU[f.FaceCell_1].rov -= info.TimeStep * flux[2]/_grid.cells[f.FaceCell_1].CellVolume;
			deltaU[f.FaceCell_1].row -= info.TimeStep * flux[3]/_grid.cells[f.FaceCell_1].CellVolume;
			deltaU[f.FaceCell_1].roE -= info.TimeStep * flux[4]/_grid.cells[f.FaceCell_1].CellVolume;*/

			if (!f.isExternal) {
				U[f.FaceCell_2].ro += info.TimeStep * flux[0]/_grid.cells[f.FaceCell_2].CellVolume;
				U[f.FaceCell_2].rou += info.TimeStep * flux[1]/_grid.cells[f.FaceCell_2].CellVolume;
				U[f.FaceCell_2].rov += info.TimeStep * flux[2]/_grid.cells[f.FaceCell_2].CellVolume;
				U[f.FaceCell_2].row += info.TimeStep * flux[3]/_grid.cells[f.FaceCell_2].CellVolume;
				U[f.FaceCell_2].roE += info.TimeStep * flux[4]/_grid.cells[f.FaceCell_2].CellVolume;

				//TO DO remove
				/*deltaU[f.FaceCell_2].ro += info.TimeStep * flux[0]/_grid.cells[f.FaceCell_2].CellVolume;
				deltaU[f.FaceCell_2].rou += info.TimeStep * flux[1]/_grid.cells[f.FaceCell_2].CellVolume;
				deltaU[f.FaceCell_2].rov += info.TimeStep * flux[2]/_grid.cells[f.FaceCell_2].CellVolume;
				deltaU[f.FaceCell_2].row += info.TimeStep * flux[3]/_grid.cells[f.FaceCell_2].CellVolume;
				deltaU[f.FaceCell_2].roE += info.TimeStep * flux[4]/_grid.cells[f.FaceCell_2].CellVolume;*/
			};

			//Add up flux to total residual
			/*double res = 0;
			for (std::set<int>::iterator it = indxs.begin(); it != indxs.end(); it++) {
				for (int k = 0; k<5; k++) res += (deltaU[*it][k] * deltaU[*it][k]);
			};
			info.Residual += sqrt(res);*/
			info.Residual = 0;
		};

		stepInfo = info;
		totalTime += stepInfo.TimeStep;
	};

	void SetInitialConditionsBlasius(double x0, double U_inf, double rho, double rhoE)
	{
		std::vector<Cell*> cells = _grid.cells.getLocalNodes();
	
		for (int i = 0; i<cells.size(); i++) {      
			Cell& cell = *cells[i];
			Vector P = cell.CellCenter;
			ConservativeVariables var;
			//at first density
			var.ro = rho;
			var.roE = rhoE;
			var.row = 0;		//if this component is zero!
			//then velocities
			if(P.x<=x0){
				var.rou = rho*U_inf;
				var.rov = 0;				
			}else{
				double y_BL = P.y*sqrt(U_inf*rho/(medium.Viscosity*(P.x - x0)));
				if(y_BL<3.5)//boundary flow
				{
					var.rou = rho*U_inf*y_BL/3.5;
					var.rov = rho*U_inf*y_BL*P.y/(7*(P.x - x0));
				}else		//outer area
				{
					var.rou = rho*U_inf;
					var.rov = 0;
				};
			};
			U[cell.GlobalIndex] = var;
		};
	
		return;
	};	
	
	//Compute viscous flux for each cell face
	void ComputeViscousFluxes() {
		//return;
		//Obtain all local faces
		std::vector<Face*> faces = _grid.faces.getLocalNodes();				

		//Compute viscous flux for each face
		double viscosity = medium.Viscosity;	
		double s_viscosity = -2.0/3.0 * viscosity;		//second viscosity
		for (int i = 0; i<faces.size(); i++) {
			Face& f = *faces[i];
			std::vector<double> vflux(5, 0);
			std::vector<Vector> grad = gradsFaces[f.GlobalIndex];
			ConservativeVariables UL = U[f.FaceCell_1]; 
            ConservativeVariables UR; 
            if (f.isExternal) { 
                UR = GetDummyCellValues(UL, f); 
            } else { 
                UR = U[f.FaceCell_2]; 
            };
			
			//Here we compute actual flux

			//required variables on face
			double rho = 0.5*(UL[0] + UR[0]);
			double u = 0.5*(UL[1] + UR[1])/rho;
			double v = 0.5*(UL[2] + UR[2])/rho;
			double w = 0.5*(UL[3] + UR[3])/rho;
			double E = 0.5*(UL[4] + UR[4])/rho;

			//gradients of required variables
			Vector du, dv, dw, dE, dT;
			du.x = (grad[1].x - u*grad[0].x)/rho;
			du.y = (grad[1].y - u*grad[0].y)/rho;
			du.z = (grad[1].z - u*grad[0].z)/rho;

			dv.x = (grad[2].x - v*grad[0].x)/rho;
			dv.y = (grad[2].y - v*grad[0].y)/rho;
			dv.z = (grad[2].z - v*grad[0].z)/rho;

			dw.x = (grad[3].x - w*grad[0].x)/rho;
			dw.y = (grad[3].y - w*grad[0].y)/rho;
			dw.z = (grad[3].z - w*grad[0].z)/rho;

			dE.x = (grad[4].x - E*grad[0].x)/rho;
			dE.y = (grad[4].y - E*grad[0].y)/rho;
			dE.z = (grad[4].z - E*grad[0].z)/rho;

			double R = 272.0;	//specific gas constant
			double K = 1.0;		//heat conduction coefficient
			dT.x = (medium.Gamma - 1)*(dE.x - u*du.x - v*dv.x - w*dw.x)/R;
			dT.y = (medium.Gamma - 1)*(dE.y - u*du.y - v*dv.y - w*dw.y)/R;
			dT.z = (medium.Gamma - 1)*(dE.z - u*du.z - v*dv.z - w*dw.z)/R;

			//stress elements
			double tau_diagonal = s_viscosity*(du.x + dv.y + dw.z);
			double tau_xx = tau_diagonal + 2*viscosity*du.x;
			double tau_yy = tau_diagonal + 2*viscosity*dv.y;
			double tau_zz = tau_diagonal + 2*viscosity*dw.z;
			double tau_xy = viscosity*(du.y + dv.x);
			double tau_xz = viscosity*(du.z + dw.x);
			double tau_yz = viscosity*(dv.z + dw.y);

			//work of viscous stresses and heat conduction
			Vector Thetta;
			Thetta.x = u*tau_xx + v*tau_xy + w*tau_xz + K*dT.x;
			Thetta.y = u*tau_xy + v*tau_yy + w*tau_yz + K*dT.y;
			Thetta.z = u*tau_xz + v*tau_yz + w*tau_zz + K*dT.z;

			//viscous fluxes
			vflux[1] = f.FaceNormal.x*tau_xx + f.FaceNormal.y*tau_xy + f.FaceNormal.z*tau_xz;			
			vflux[2] = f.FaceNormal.x*tau_xy + f.FaceNormal.y*tau_yy + f.FaceNormal.z*tau_yz;			
			vflux[3] = f.FaceNormal.x*tau_xz + f.FaceNormal.y*tau_yz + f.FaceNormal.z*tau_zz;			
			vflux[4] = f.FaceNormal*Thetta;
			
			//
			for (int k = 0; k<5; k++) vflux[k] *= f.FaceSquare;

			//Add viscous fluxes			
			fluxes[f.GlobalIndex] += vflux; //Тут по идее он должен быть умноже на площадь грани так что обрати внимание на свои формулы
		};
	};


	//Compute gradient at cell given its neighbours
	std::vector<Vector> ComputeGradInCell(Cell& cell, std::set<int> nIdx) {
		//Prepare data containers
		std::vector<Vector> grad(5);		
		std::vector<int> neighbours;
		for (std::set<int>::iterator it = nIdx.begin(); it!=nIdx.end(); it++) {
			neighbours.push_back(*it);			
		};

		//Build matrix (depends only on grid)
		alglib::real_2d_array matrix;
		matrix.setlength(nIdx.size(), 3);
		for (int i = 0; i<neighbours.size(); i++) {
			Cell& c = _grid.cells[neighbours[i]];
			Vector dr = c.CellCenter - cell.CellCenter;
			matrix[i][0] = dr.x;
			matrix[i][1] = dr.y;
			matrix[i][2] = dr.z;
		};

		//Solve for gradient with chosen right hand side
		alglib::real_1d_array rhs;
		rhs.setlength(neighbours.size());
		alglib::real_1d_array x;		
		x.setlength(3);
		alglib::ae_int_t info;
		alglib::lsfitreport rep;
		      
		for (int k = 0; k<5; k++) {			
			//Construct right hand side
			double uCell = U[cell.GlobalIndex][k];
			for (int i = 0; i<neighbours.size(); i++) {
				Cell& c = _grid.cells[neighbours[i]];
				double dU = U[c.GlobalIndex][k] - uCell;
				rhs[i] = dU;
			};

			//Solve system for x in linear least sense
			// Linear fitting without weights
			alglib::lsfitlinear(rhs, matrix, info, x, rep);
			grad[k] = Vector(x[0], x[1], x[2]);
		};

		return grad;
	};

	//Function to compute gradients at every cell
	void ComputeGradients() {
		gradsCells.clear();
		std::vector<Cell*> cells = _grid.cells.getLocalNodes();

		//First we build kd-Tree for every cell center		
		alglib::ae_int_t nx = 3;
		alglib::ae_int_t ny = 1;
		alglib::ae_int_t normtype = 2;		
		alglib::real_2d_array a;  		
		a.setlength(cells.size(), nx + ny);				
		
		for (int i = 0; i<cells.size(); i++) {						
			Cell& cell = *cells[i];

			//Fill in kdtree information
			a[i][0] = cell.CellCenter.x;
			a[i][1] = cell.CellCenter.y;
			a[i][2] = cell.CellCenter.z;
			a[i][3] = cell.GlobalIndex;			
		};
	
		// Build tree	
		alglib::kdtreebuild(a, nx, ny, normtype, kdtCells);		

		//Now for each local cell obtain nearest neighbour cells
		int kNeighbours = _grid.gridInfo.GridDimensions;		
		alglib::real_1d_array x; 
		x.setlength(nx);
		alglib::real_2d_array xy; 		
		xy.setlength(kNeighbours, nx+ny);
		gradsCells.clear();		
		cellNeigbours.clear();
		for (int i = 0; i<cells.size(); i++) {						
			Cell& cell = *cells[i];

			//Get nearest neighbours
			x[0] = cell.CellCenter.x;
			x[1] = cell.CellCenter.y;
			x[2] = cell.CellCenter.z;			
			alglib::ae_int_t isFound = alglib::kdtreequeryknn(kdtCells, x, kNeighbours);
			alglib::kdtreequeryresultsxy(kdtCells, xy);
			for (int k = 0; k<kNeighbours; k++) cellNeigbours[cell.GlobalIndex].insert((int)xy[k][3]);			
		};

		////Using least squares compute gradient at cell center
		for (int i = 0; i<cells.size(); i++) {			
			Cell& cell = *cells[i];			
			gradsCells[cell.GlobalIndex] = ComputeGradInCell(cell, cellNeigbours[cell.GlobalIndex]);
		};		
						
		//Compute gradients at face centers
		//Obtain all local faces
		std::vector<Face*> faces = _grid.faces.getLocalNodes();		

		for (int i = 0; i<faces.size(); i++) {
			Face& face = *faces[i];
			//For external faces just get gradient in nearest cell
			if (face.isExternal) {
				gradsFaces[face.GlobalIndex] = gradsCells[face.FaceCell_1];			
			} else {
				gradsFaces[face.GlobalIndex] = (gradsCells[face.FaceCell_1] + gradsCells[face.FaceCell_2]);			
				gradsFaces[face.GlobalIndex] /= 2.0;
			};
		};
	};

	//Compute wall distances for every cell
	void ComputeWallDistances() {		
		//If no walls defined
		if (_wallFaces.size() == 0) {
			std::cout<<"No walls defined. Distance computation is omitted.\n";
			return;
		};		

		// Build alglib kdtree for all wall face centers
		alglib::ae_int_t nx = 3;
		alglib::ae_int_t ny = 1;
		alglib::ae_int_t normtype = 2;		
		alglib::real_2d_array a;  		
		a.setlength(_wallFaces.size(), nx + ny);		

		int i = 0;
		for (std::map<int, FaceWallInfo>::iterator it = _wallFaces.begin(); it != _wallFaces.end(); it++) {			
			// Clear all wall info for faces
			it->second.layerCells.clear();

			//Fill in kdtree information
			a[i][0] = _grid.faces[it->first].FaceCenter.x;
			a[i][1] = _grid.faces[it->first].FaceCenter.y;
			a[i][2] = _grid.faces[it->first].FaceCenter.z;
			a[i][3] = _grid.faces[it->first].GlobalIndex;
			i++;
		};
	
		// Build tree	
		alglib::kdtreebuild(a, nx, ny, normtype, kdtWall);		

		//Now for each local cell obtain distance to closest wall face and other properties
		std::vector<Cell*> cells = _grid.cells.getLocalNodes();
		alglib::real_1d_array x; 
		x.setlength(nx);
		alglib::real_2d_array xy; 		
		xy.setlength(1, nx+ny);
		for (int i = 0; i < cells.size(); i++) {			
			x[0] = cells[i]->CellCenter.x;
			x[1] = cells[i]->CellCenter.y;
			x[2] = cells[i]->CellCenter.z;			
			alglib::ae_int_t isFound = alglib::kdtreequeryknn(kdtWall, x, 1);
			alglib::kdtreequeryresultsxy(kdtWall, xy);						
			Vector radiusVector = Vector(xy[0][0], xy[0][1], xy[0][2]);			
			radiusVector -= cells[i]->CellCenter;
			int faceIndex = xy[0][3];
			_wallInfo[cells[i]->GlobalIndex] = CellWallInfo();		
			_wallInfo[cells[i]->GlobalIndex].wallFaceIndex = faceIndex;
			_wallInfo[cells[i]->GlobalIndex].distance = radiusVector.mod();
			_wallInfo[cells[i]->GlobalIndex].angle = radiusVector.mod();
			_wallFaces[faceIndex].layerCells.push_back(cells[i]->GlobalIndex);
		};
	};
	

	void ComputeDisplacementThickness(int wallFaceInd, double U_inf)
	{
		//layer cells must be sorted in ascending order by distance
		std::vector<int> CellSequence = _wallFaces[wallFaceInd].layerCells;
		
		//write required function
		function f(CellSequense.size()+1);
		f.x[0] = 0;		//values on wall
		f.value[0] = 0;
		for(int i=0; i<CellSequense.size(); i++)
		{
			f.x[i+1] = _wallInfo[CellSequence[i]].distance;
			f.value[i+1] = 1 - U[CellSequence[i]].rou/(U[CellSequence[i]].ro*U_inf);
		};
		//compute integral by trapezium formula
		_wallFaces[wallFaceInd].boundaryLayerThickness = f.TrapeziumIntegral()
	};

	//Boundary conditions main processing function
	ConservativeVariables GetDummyCellValues(ConservativeVariables U, const Face& face) {
		//Determine boundary condition type
		BoundaryCondition* bc = _boundaryConditions[face.BCMarker];
		return bc->getDummyValues(U, face); 
	};

	//


	//Conversion functions and postprocessing
	ConservativeVariables PrimitiveToConservativeVariables(Vector v, double p, double T, MediumProperties medium) {
		ConservativeVariables res;
		double e = medium.Cv * T;
		double v2 = v*v; 
		res.ro = p / (e * (medium.Gamma - 1));
		res.rou = v.x * res.ro;
		res.rov = v.y * res.ro;
		res.row = v.z * res.ro;
		res.roE = (e + v2 / 2.0)* res.ro;
		return res;
	};

	Vector GetVelocity(const ConservativeVariables& celldata) {		
		double ro = celldata.ro;
		double vx = celldata.rou/celldata.ro;
		double vy = celldata.rov/celldata.ro;
		double vz = celldata.row/celldata.ro;
		return Vector(vx, vy, vz);
	};

	double GetPressure(const ConservativeVariables& celldata) {
		double ro = celldata.ro;
		double vx = celldata.rou/celldata.ro;
		double vy = celldata.rov/celldata.ro;
		double vz = celldata.row/celldata.ro;
		double E = celldata.roE/celldata.ro;
		double P = (medium.Gamma-1.0) * ro * (E - (vx*vx+vy*vy+vz*vz)/2.0);
		return P;
	};

	double GetTemperature(const ConservativeVariables& celldata) {
		double ro = celldata.ro;
		double vx = celldata.rou/celldata.ro;
		double vy = celldata.rov/celldata.ro;
		double vz = celldata.row/celldata.ro;
		double E = celldata.roE/celldata.ro;
		double T = (E - (vx*vx+vy*vy+vz*vz)/2.0) / (medium.Cv);
		return T;
	};

	double GetSoundSpeed(const ConservativeVariables& celldata) {
		double ro = celldata.ro;
		double vx = celldata.rou/celldata.ro;
		double vy = celldata.rov/celldata.ro;
		double vz = celldata.row/celldata.ro;
		double E = celldata.roE/celldata.ro;
		double P = (medium.Gamma-1.0) * ro * (E - (vx*vx+vy*vy+vz*vz)/2.0);
		double c = sqrt(medium.Gamma * P / ro);
		return c;
	};


	//Save current model state to techplot file
	void SaveToTechPlot(std::string fname) {
		std::ofstream ofs(fname);
		//TO DO unify		
		//Header
		ofs<<"VARIABLES= \"X\", \"Y\", \"Rho\", \"u\", \"v\", \"w\", \"E\", \"T\", \"P\", \"distance\"\n";
		ofs<<"ZONE T=\"D\"\n";
		ofs<<"N=" << _grid.nodes.size() << ", E=" << _grid.cells.size() <<", F=FEBLOCK, ET=QUADRILATERAL\n";
		ofs<<"VARLOCATION = (NODAL, NODAL, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED, CELLCENTERED)\n";

		//Map all node global indexes to natural numbers
		std::map<int,int> toNaturalIndex;
		std::set<int> nodeIndexes = _grid.nodes.getAllIndexes();
		int counter = 1;
		for (std::set<int>::iterator it = nodeIndexes.begin(); it != nodeIndexes.end(); it++) toNaturalIndex[*it] = counter++;


		//Access local nodes, faces, cells and flow data
		std::vector<Node*> nodes = _grid.nodes.getLocalNodes();
		std::vector<Face*> faces = _grid.faces.getLocalNodes();
		std::vector<Cell*> cells = _grid.cells.getLocalNodes();
		std::vector<ConservativeVariables*> data = U.getLocalNodes();

		//Nodes coordinates
		//X
		for (int i = 0; i<nodes.size(); i++) {
			ofs<<nodes[i]->P.x<<"\n";
		};

		//Y
		for (int i = 0; i<nodes.size(); i++) {
			ofs<<nodes[i]->P.y<<"\n";
		};

		////Solution data
		////Rho
		for (int i = 0; i<cells.size(); i++) {
			int idx = cells[i]->GlobalIndex;
			ofs<<U[idx].ro<<"\n";
		};
				
		//u
		for (int i = 0; i<cells.size(); i++) {
			int idx = cells[i]->GlobalIndex;
			ofs<<U[idx].rou/U[idx].ro<<"\n";
		};
		//v
		for (int i = 0; i<cells.size(); i++) {
			int idx = cells[i]->GlobalIndex;
			ofs<<U[idx].rov/U[idx].ro<<"\n";		
		};
		//w
		for (int i = 0; i<cells.size(); i++) {
			int idx = cells[i]->GlobalIndex;
			ofs<<U[idx].row/U[idx].ro<<"\n";
		};
		//E
		for (int i = 0; i<cells.size(); i++) {
			int idx = cells[i]->GlobalIndex;
			ofs<<U[idx].roE/U[idx].ro<<"\n";
		};
		//T
		for (int i = 0; i<cells.size(); i++) {
			int idx = cells[i]->GlobalIndex;
			double ro = U[idx].ro;
			double vx = U[idx].rou/U[idx].ro;
			double vy = U[idx].rov/U[idx].ro;
			double vz = U[idx].row/U[idx].ro;
			double E = U[idx].roE/U[idx].ro;
			double T = (E - (vx*vx+vy*vy+vz*vz)/2.0) / (medium.Cv);
			ofs<<T<<"\n";
		};
		//P
		for (int i = 0; i<cells.size(); i++) {
			int idx = cells[i]->GlobalIndex;
			double ro = U[idx].ro;
			double vx = U[idx].rou/U[idx].ro;
			double vy = U[idx].rov/U[idx].ro;
			double vz = U[idx].row/U[idx].ro;
			double E = U[idx].roE/U[idx].ro;
			double P = (medium.Gamma-1.0) * ro * (E - (vx*vx+vy*vy+vz*vz)/2.0);
			ofs<<P<<"\n";
		};		

		//distance
		for (int i = 0; i<cells.size(); i++) {
			int idx = cells[i]->GlobalIndex;
			if (_wallInfo.size() != 0) {
				ofs<<_wallInfo[idx].distance<<"\n";
			} else {
				ofs<<0<<"\n";
			};
		};

		//Connectivity list for each cell
		for (int i = 0; i<cells.size(); i++) {
			ofs<<toNaturalIndex[cells[i]->Nodes[0]]<<" "
				<<toNaturalIndex[cells[i]->Nodes[1]]<<" "
				<<toNaturalIndex[cells[i]->Nodes[2]]<<" "
				<<toNaturalIndex[cells[i]->Nodes[3]]<<"\n";
		};

		ofs.close();
		return;

	};


	//Save solution in a file
	void SaveSolution(std::string fname = "SavedSolution.txt")
	{
		std::ofstream ofs(fname);
		std::vector<Cell*> cells = _grid.cells.getLocalNodes();
		ofs << totalTime << "\n";			
		for(int i=0; i<cells.size(); i++)
		{
			int Index = cells[i]->GlobalIndex;	//Global Index of cell
			ofs << Index << " " << U[Index].ro << " " << U[Index].roE << " " << U[Index].rou << " ";
			ofs << U[Index].rov << " " << U[Index].row << "\n";
		};
		ofs.close();
	};
	
	//Load a solution file
	void LoadSolution(std::string fname = "SavedSolution.txt")
	{
		std::ifstream ifs(fname);
		ifs >> totalTime;				
		while(!ifs.eof())
		{
			int GlobalInd;
			ifs >> GlobalInd;
			ifs >> U[GlobalInd].ro >> U[GlobalInd].roE >> U[GlobalInd].rou >> U[GlobalInd].rov >> U[GlobalInd].row;
		};
		ifs.close();
	};

};


#endif